/**
 * to count the matrix, if use this tool before GSNP, 
 * then you can turn off the count_matrix generation step
 * */

#include "soap_snp.h"


void generateCountMatrix(ubit64_t count_matrix[COUNT_MATRIX_SIZE], 
			std::ifstream &alignment, Genome* genome) {
	// Read Alignment files
	Soap_format soap;
	map<Chr_name, Chr_info*>::iterator current_chr;
	current_chr = genome->chromosomes.end();
	ubit64_t ref(0);
	std::string::size_type coord;

	unsigned int count = 0;

	for(std::string line; getline(alignment, line);) {
		std::istringstream s(line);
		if( s>> soap ) {

			g_numTotalSite++;

			if(soap.position < 0) {
				continue;
			}
			//cerr<<soap<<endl;
			// In the overloaded "+" above, soap.position will be substracted by 1 so that coordiates start from 0
			if (current_chr == genome->chromosomes.end() || current_chr->first != soap.chr_name) {
				current_chr = genome->chromosomes.find(soap.chr_name);
				cout << "here 1..." << endl;
				if(current_chr == genome->chromosomes.end()) {
					for(map<Chr_name, Chr_info*>::iterator test = genome->chromosomes.begin();test != genome->chromosomes.end();test++) {
						cerr<<'!'<<(test->first)<<'!'<<endl;
					}
					cerr<<"Assertion Failed: Chromosome: !"<<soap.chr_name<<"! NOT found"<<endl;
					exit(255);
				}
			}
			else {
				;
			}
			if (soap.position + soap.read_len>=current_chr->second->length()) {
				continue;
			}
			if (soap.hit == 1) {
				for(coord = 0; coord != soap.read_len; coord++) {
					if (soap.read[coord] == 'N') {
						;
					}
					else {
						if(! (soap.position + coord<current_chr->second->length())) {
							cerr<<soap<<endl;
							cerr<<"The program found the above read has exceed the reference length:\n";
							cerr<<"The read is aligned to postion: "<<soap.position<<" with read length: "<<soap.read_len<<endl;
							cerr<<"Reference: "<<current_chr->first<<" FASTA Length: "<<current_chr->second->length()<<endl;
							exit(255);
						}
						ref = current_chr->second->get_bin_base(soap.position+coord);
						if ( (ref&12) !=0 ) {
							// This is an N on reference or a dbSNP which should be excluded from calibration
							;
						}
						else {

							int idx = 0;
							if(soap.strand=='+') {
								// forward strand
								count_matrix[(((ubit64_t)soap.qual[coord])<<12) | (coord<<4) | ((ref&0x3)<<2) | ((soap.read[coord]>>1)&3)] += 1;
								idx = (((ubit64_t)soap.qual[coord])<<12) | (coord<<4) | ((ref&0x3)<<2) | ((soap.read[coord]>>1)&3);
							}
							else {
								// reverse strand
								count_matrix[(((ubit64_t)soap.qual[coord])<<12) | ((soap.read_len-1-coord)<<4) | ((ref&0x3)<<2) | ((soap.read[coord]>>1)&3)] += 1;
								idx = (((ubit64_t)soap.qual[coord])<<12) | ((soap.read_len-1-coord)<<4) | ((ref&0x3)<<2) | ((soap.read[coord]>>1)&3);
							}
						}
					}
				}
			}
		}
	}
}

int main(int argc, char** argv) {
	/*parameter input*/
	cout << "This is count_matrx" << endl;
	if(argc != 3 && argc != 4) {
		cout << "Please specific three parameters:" << endl;
		cout << "1. The reference sequence file." << endl;
		cout << "2. The SOAP alignment file." << endl;
		cout << "3. The dbSNP file." << endl;
		exit(-1);
	}
	const char* refFileName = argv[1];
	const char* alignmentFileName = argv[2];
	cout << "Reference file: " << refFileName << endl;
	cout << "Alignment file: " << alignmentFileName << endl;
	bool hasDbsnp = false;
	char* dbsnpFileName;
	if(argc == 4) {
		dbsnpFileName = argv[3];
		hasDbsnp = true;
		cout << "dbSNP file: " << dbsnpFileName << endl;
	}

	/*parameter initialization and open files*/
	ifstream refFile(refFileName, ifstream::in);
	ifstream alignmentFile(alignmentFileName, ifstream::in);
	ifstream snpFile;
	if(hasDbsnp) {
		snpFile.open(dbsnpFileName);
	}
	Genome* genome = new Genome(refFile, snpFile);

	/*generate the counting matrix*/
	ubit64_t count_matrix[COUNT_MATRIX_SIZE] = {0};
	generateCountMatrix(count_matrix, alignmentFile, genome);


	/*output the couting matrix as the result*/
	cout << "------------------------" << endl;
	string outFileName = string(alignmentFileName) + ".cmat";
	cout << "Output file: " << outFileName << endl;
	FILE* outFile;
	safe_fopen(&outFile, outFileName.c_str(), "wb");
	safe_fwrite(count_matrix, sizeof(ubit64_t), COUNT_MATRIX_SIZE, outFile);
	safe_fclose(outFile);
	cout << "Done." << endl;

	/*cleanup*/
	refFile.close();
	alignmentFile.close();

	return EXIT_SUCCESS;
}

